##############################################################################################################################
#### Title: CPS Replication Files for “Legacy of Resisting State Repression”
#### Author: Sangyong Son
##############################################################################################################################

##############################################################################################################################
#### Load Packages 
##############################################################################################################################
library(readxl); library(ggplot2); library(dplyr); library(stargazer); library(jtools); library(broom); library(ggthemes);
library(gridExtra); library(ggplot2); library(ggeffects); library(tibble); library(broom); library(lmtest);
library(gridtext); library(grid); library(cobalt); library(MatchIt); library(estimatr)

##############################################################################################################################
#### Load Data
##############################################################################################################################
Data <- readRDS("~/..../Son_Data.rds")


##############################################################################################################################
#### Table 3: Selection into responses to repression and their correlation with post-repression political activism
##############################################################################################################################
Data <- Data %>%
  mutate(
    resistance1 = case_when(
      Collective.Action == 1 ~ 0,
      Collective.Action %in% c(2, 3) ~ 1,
      TRUE ~ NA_real_
    ),
    resistance2 = case_when(
      Collective.Action == 2 ~ 0,
      Collective.Action == 3 ~ 1,
      TRUE ~ NA_real_   
    ),
    
    Collective.Action = factor(
      Collective.Action,
      levels = c(1, 2, 3),
      labels = c("Non-resisters", "Individualized resisters", "Institutionalized resisters")
    )
  )

## Resistance1: Participation in Resistance (Resisters vs. Nonresisters)
fit1 <- lm(resistance1~Age+Sex+Marriage+Employment+Economic.Status+Location.2+Education+University.Student+
             CA1+CA6+CA7+CA8+CA9, data=Data); summary(fit1)

## Resistance2: Participation in Institutionalized Resistance (Individualized Resisters vs. Institutionalized Resisters)
fit2 <- lm(resistance2~Age+Sex+Marriage+Employment+Economic.Status+Location.2+Education+University.Student+
             CA1+CA6+CA7+CA8+CA9, data=Data); summary(fit2)

fit3 <- lm(Political.Participation.Index~Collective.Action+Age+Sex+Marriage+Employment+Economic.Status+Location.2+
             Education+University.Student+CA1+CA6+CA7+CA8+CA9,data = Data); summary(fit3)

fit4 <- lm(Contribution_A~Collective.Action+Age+Sex+Marriage+Employment+Economic.Status+Location.2+Education+University.Student+
             +CA1+CA6+CA7+CA8+CA9, data = Data); summary(fit4)

fit5 <- lm(Contribution_B~Collective.Action+Age+Sex+Marriage+Employment+Economic.Status+Location.2+Education+University.Student+
             +CA1+CA6+CA7+CA8+CA9, data = Data); summary(fit5)

## Table 3 in the manuscript
stargazer(fit1, fit2, fit3, fit4, fit5)  

##############################################################################################################################
#### Matching  
##############################################################################################################################
Data <- readRDS("~/..../Son_Data.rds")

##############################################################################################################################
#### (1) Effects of Participation in Resistance: Resisters vs. Nonresisters 
##############################################################################################################################
Data1 <- Data
Data1$Collective.Action <- car::recode(Data1$Collective.Action, "1=0; c(2,3)=1")

Data1 <- Data1 %>% dplyr::select(
  Age, Sex, Education, Marriage, Employment, Location.2, CA1, CA6, CA7, CA8, CA9, CA10,
  University.Student, Economic.Status, Collective.Action,
  Contribution_A, Contribution_B, Political.Participation.Index,
  Low.Participation.Index, High.Participation.Index, Career, Collective.Action.Violence
)

resistance.subclass <- matchit(
  Collective.Action ~ Age + Sex + Employment + Economic.Status + Marriage +
    Location.2 + Education + University.Student,
  data = Data1, method = "subclass"
)

new.names <- c(
  Economic.Status = "Economic status",
  Location.2 = "Location",
  University.Student = "University student",
  CA1 = "Night school movement",
  CA6 = "Student movement",
  CA7 = "Christian association",
  CA8 = "Labor movement",
  CA9 = "Women's movement"
)

b1 <- love.plot(
  Collective.Action ~ Age + Sex + Employment + Marriage + Location.2 +
    Education + University.Student + CA1 + CA6 + CA7 + CA8 + CA9,
  data = Data1, estimand = "ATT",
  weights = list(w1 = get.w(resistance.subclass)),
  var.order = "unadjusted", abs = TRUE, line = TRUE,
  thresholds = c(m = .1), var.names = new.names, limits = c(0, .82),
  colors = c("red","blue"), shapes = c("triangle filled","circle filled"),
  sample.names = c("Unadjusted","Subclass")
) +
  labs(title = "Effects of joining resistance \n Resisters vs non-resisters") +
  theme(legend.position = c(.75, .13),
        legend.box.background = element_rect(),
        legend.box.margin = margin(1, 1, 1, 1),
        plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
        axis.title.x = element_text(size = 15, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        axis.text    = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.text  = element_text(size = 15))


m.data1 <- match.data(resistance.subclass)
m.data1$Collective.Action <- car::recode(m.data1$Collective.Action, "0='Non-resisters'; 1='Resisters'")
m.data1$Collective.Action <- as.factor(m.data1$Collective.Action)
m.data1$Collective.Action <- relevel(m.data1$Collective.Action, "Non-resisters")

fit6 <- lm(Political.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
             Employment + Economic.Status + Location.2 + Education + University.Student +
             CA1 + CA6 + CA7 + CA8 + CA9,
           data = m.data1, weights = weights)

fit7 <- glm(Contribution_A ~ Collective.Action + Age + Sex + Marriage +
              Employment + Economic.Status + Location.2 + Education + University.Student,
            data = m.data1, weights = weights, family = binomial())

fit8 <- glm(Contribution_B ~ Collective.Action + Age + Sex + Marriage +
              Employment + Economic.Status + Location.2 + Education + University.Student,
            data = m.data1, weights = weights, family = binomial())

fit7.ctrls <- lm(Contribution_A ~ Collective.Action + Age + Sex + Marriage +
             Employment + Economic.Status + Location.2 + Education + University.Student+
             CA1 + CA6 + CA7 + CA8 + CA9,
           data = m.data1, weights = weights)

fit8.ctrls <- lm(Contribution_B ~ Collective.Action + Age + Sex + Marriage +
             Employment + Economic.Status + Location.2 + Education + University.Student+
             CA1 + CA6 + CA7 + CA8 + CA9,
           data = m.data1, weights = weights)


##############################################################################################################################
#### (2) Effects of Joining Institutionalized Resistance: Individualized vs. Institutionalized Resisters   
##############################################################################################################################
Data2 <- Data
Data2 <- dplyr::filter(Data2, Collective.Action %in% c(2, 3))
Data2$Collective.Action <- car::recode(Data2$Collective.Action, "3=1; 2=0")

Data2 <- Data2 %>% dplyr::select(
  Age, Sex, Education, Marriage, Employment, Location.2, CA1, CA6, CA7, CA8, CA9, CA10,
  University.Student, Economic.Status, Collective.Action,
  Contribution_A, Contribution_B, Political.Participation.Index,
  Low.Participation.Index, High.Participation.Index, Career, Collective.Action.Violence
)

institutionalized.resistance.subclass <- matchit(
  Collective.Action ~ Age + Sex + Marriage + Employment + Economic.Status +
    Location.2 + Education + University.Student+
    +CA1+CA6+CA7+CA8+CA9,
  data = Data2, method = "subclass"
)

new.names <- c(
  Economic.Status = "Economic status",
  Location.2 = "Location",
  University.Student = "University student",
  CA1 = "Night school movement",
  CA6 = "Student movement",
  CA7 = "Christian association",
  CA8 = "Labor movement",
  CA9 = "Women's movement"
)

b2 <- love.plot(
  Collective.Action ~ Age + Sex + Employment + Economic.Status + Marriage +
    Location.2 + Education + University.Student + CA1 + CA6 + CA7 + CA8 + CA9,
  data = Data2, estimand = "ATT",
  weights = list(w1 = get.w(institutionalized.resistance.subclass)),
  var.order = "unadjusted",
  abs = TRUE,
  line = TRUE,
  thresholds = c(m = .1),
  var.names = new.names,
  colors = c("red","blue"),
  shapes = c("triangle filled","circle filled"),
  sample.names = c("Unadjusted","Subclass"),
  limits = c(0, .82)
) +
  labs(title = "Effects of joining institutionalized resistance \n Individualized vs institutionalized resisters") +
  theme(legend.position = c(.75, .13),
        legend.box.background = element_rect(),
        legend.box.margin = margin(1, 1, 1, 1),
        plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
        axis.title.x = element_text(size = 15, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        axis.text    = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.text  = element_text(size = 15))

grid.arrange(b1, b2, ncol = 2) ## Figure 3 in the Supplementary Information 

m.data2 <- match.data(institutionalized.resistance.subclass)
m.data2$Collective.Action <- car::recode(m.data2$Collective.Action, "0='Individual'; 1='Institutional'")
m.data2$Collective.Action <- as.factor(m.data2$Collective.Action)

fit9  <- lm(Political.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
              Employment + Economic.Status + Location.2 + Education + University.Student +
              CA1 + CA6 + CA7 + CA8 + CA9,
            data = m.data2, weights = weights)

fit10 <- glm(Contribution_A ~ Collective.Action + Age + Sex + Marriage +
               Employment + Economic.Status + Location.2 + Education + University.Student,
             data = m.data2, weights = weights, family = binomial())

fit11 <- glm(Contribution_B ~ Collective.Action + Age + Sex + Marriage +
               Employment + Economic.Status + Location.2 + Education + University.Student,
             data = m.data2, weights = weights, family = binomial())

fit10.ctrls <- lm(Contribution_A ~ Collective.Action + Age + Sex + Marriage +
               Employment + Economic.Status + Location.2 + Education + University.Student+
                 CA1 + CA6 + CA7 + CA8 + CA9,
             data = m.data2, weights = weights)

fit11.ctrls <- lm(Contribution_B ~ Collective.Action + Age + Sex + Marriage +
               Employment + Economic.Status + Location.2 + Education + University.Student+
                 CA1 + CA6 + CA7 + CA8 + CA9,
             data = m.data2, weights = weights)

#### Table 3 in the Supplementary Information 
stargazer(fit6, fit7.ctrls, fit8.ctrls, fit9, fit10.ctrls, fit11.ctrls)

#######################################################################################################################
###### Figure 2 
#######################################################################################################################
dat1  <- ggpredict(fit6, terms = "Collective.Action")
dat11 <- ggpredict(fit9, terms = "Collective.Action")

pp_theme <- theme_classic() +
  theme(
    plot.title  = element_text(hjust = 0.5, size = 18, face = "bold"),
    axis.title.x= element_text(size = 18, face = "bold", vjust = -0.6),
    axis.title.y= element_text(size = 15, face = "bold"),
    axis.text.x = element_text(size = 16, face = "bold"),
    axis.text.y = element_text(size = 14)
  )

p1 <- plot(dat1, ci.style = "errorbar", dot.size = 4) +
  pp_theme +
  ggtitle("Predicted values of \n political participation (0-8)") +
  xlab("") + ylab("") +
  scale_x_discrete(limits = c("Non-resisters","Resisters"))

p2 <- plot(dat11, ci.style = "errorbar", dot.size = 4) +
  pp_theme +
  ggtitle("Predicted values of \n political participation (0-8)") +
  xlab("") + ylab("") +
  scale_x_discrete(limits = c("Individualized resisters","Institutionalized resisters"))

# Figure 2 in the manuscript 
grid.arrange(p1, p2, ncol = 2)  


#######################################################################################################################
###### Figure 3 
#######################################################################################################################
## When visualizing results for public good games, I use GLM models without social network control variables due to separation issues. 
## Specifically, I use fit7, fit8, fit10, and fit11.

pp_theme <- theme_classic() +
  theme(
    plot.title  = element_text(hjust = 0.5, size = 18, face = "bold"),
    axis.title.x= element_text(size = 18, face = "bold", vjust = -0.6),
    axis.title.y= element_text(size = 15, face = "bold"),
    axis.text.x = element_text(size = 16, face = "bold"),
    axis.text.y = element_text(size = 14)
  )

make_plot <- function(dat, title_text, limits_vec = NULL) {
  p <- plot(dat, ci.style = "errorbar", dot.size = 4) +
    pp_theme +
    ggtitle(title_text) +
    xlab("") + ylab("")
  if (!is.null(limits_vec)) {
    p <- p + scale_x_discrete(limits = limits_vec)
  }
  p
}

#### 50% condition -------------------------------------------------------------
dat2 <- ggpredict(fit7,  terms = "Collective.Action")  ## 39/76
dat3 <- ggpredict(fit10, terms = "Collective.Action")  ## 77/98

Plot2 <- make_plot(dat2, "Effects of joining resistance",
                   limits_vec = levels(dat2$x)) +
  annotate("label", x = 2, y = 0.73, label = "76%", size = 6, fontface = "bold") +
  annotate("label", x = 1, y = 0.36, label = "39%", size = 6, fontface = "bold")

Plot3 <- make_plot(dat3, "Effects of joining institutionalized resistance",
                   limits_vec = levels(dat3$x)) +
  annotate("label", x = 2, y = 0.96, label = "98%", size = 6, fontface = "bold") +
  annotate("label", x = 1, y = 0.75, label = "77%", size = 6, fontface = "bold")

#### Figure 3 of the Manuscript 
k1 <- grid.arrange(
  Plot2, Plot3, ncol = 2, padding = unit(1.5, "line"),
  top = textGrob(
    expression(bold("Predicted values of contributors for public goods in 50% threshold condition")),
    gp = gpar(fontsize = 23)
  )
)

#### 75% condition -------------------------------------------------------------
dat5 <- ggpredict(fit8,  terms = "Collective.Action")  ## 34/69
dat6 <- ggpredict(fit11, terms = "Collective.Action")  ## 69/97

Plot5 <- make_plot(dat5, "Effects of joining resistance",
                   limits_vec = levels(dat5$x)) +
  annotate("label", x = 2, y = 0.66, label = "69%", size = 6, fontface = "bold") +
  annotate("label", x = 1, y = 0.30, label = "34%", size = 6, fontface = "bold")

Plot6 <- make_plot(dat6, "Effects of joining institutionalized resistance",
                   limits_vec = levels(dat6$x)) +
  annotate("label", x = 2, y = 0.938, label = "97%", size = 6, fontface = "bold") +
  annotate("label", x = 1, y = 0.666, label = "69%", size = 6, fontface = "bold")

#### Figure 3 of the Manuscript 
k2 <- grid.arrange(
  Plot5, Plot6, ncol = 2, padding = unit(1.5, "line"),
  top = textGrob(
    expression(bold("Predicted values of contributors for public goods in 75% threshold condition")),
    gp = gpar(fontsize = 23)  
  )
)

#######################################################################################################################
###### Table 5: Career Trajectory Analysis 
#######################################################################################################################
fit.career1 <- lm(Career ~ Collective.Action+Age+Sex+Marriage+
                    Employment+Economic.Status+Location.2+Education+University.Student+CA1+CA6+CA7+CA8+CA9, 
                  data = m.data1, weights = weights); summary(fit.career1)

fit.career2 <- lm(Career ~ Collective.Action+Age+Sex+Marriage+
                    Employment+Economic.Status+Location.2+Education+University.Student+CA1+CA6+CA7+CA8+CA9, 
                  data = m.data2, weights = weights); summary(fit.career2)

#### Table 5 of the Manuscript 
stargazer(fit.career1, fit.career2)

#######################################################################################################################
###### Robustness Checks (1)
#######################################################################################################################
controls_all <- c("Age","Sex","Marriage","Employment","Economic.Status",
                  "Location.2","Education","University.Student","CA1","CA6","CA7","CA8","CA9")

f_lin <- function(y, controls = controls_all) {
  as.formula(paste(y, "~ as.factor(Collective.Action) +", paste(controls, collapse = " + ")))
}

run_three_lm <- function(df, wvar = "weights") {
  df$.__w__ <- df[[wvar]] 
  list(
    lm(f_lin("Political.Participation.Index"), data = df, weights = .__w__),
    lm(f_lin("Contribution_A"),                data = df, weights = .__w__),
    lm(f_lin("Contribution_B"),                data = df, weights = .__w__)
  )
}

drop_students <- function(df) dplyr::filter(df, University.Student != 1)
drop_CA_any1  <- function(df) {
  dplyr::filter(df, if_all(dplyr::all_of(c("CA1","CA6","CA7","CA8","CA9")), ~ . != 1))
}


show_means <- function(df) {
  round(c(Mean_PP = mean(df$Political.Participation.Index),
          Mean_A  = mean(df$Contribution_A),
          Mean_B  = mean(df$Contribution_B)), 2)
}

############################################################################################################
# A: Joining resistance
############################################################################################################
Data1 <- m.data1
table(Data1$University.Student)    # 92

# University students drop
Data1 <- drop_students(Data1)
table(Data1$University.Student)    # 0
nrow(Data1)                        # 608

fits_A_students <- run_three_lm(Data1)
fit1r1 <- fits_A_students[[1]]; pro1r1 <- fits_A_students[[2]]; pro2r1 <- fits_A_students[[3]]

show_means(Data1)                  # 4.79, 0.52, 0.46

# Resisters (movement dummies = 1) drop
Data2 <- m.data3
nrow(Data2)                        # 700

Data2 <- drop_CA_any1(Data2)
nrow(Data2)                        # 645

fits_A_resisters <- run_three_lm(Data2)
fit2r1 <- fits_A_resisters[[1]]; pro3r1 <- fits_A_resisters[[2]]; pro4r1 <- fits_A_resisters[[3]]

show_means(Data2)                  # 4.69, 0.51, 0.45

stargazer::stargazer(fit1r1, pro1r1, pro2r1) ## Table 4 of the Supplementary Information 
stargazer::stargazer(fit2r1, pro3r1, pro4r1) ## Table 5 of the Supplementary Information

############################################################################################################
# B: Joining resistance institution
############################################################################################################
Data3 <- m.data2
nrow(Data3)                        # 318

# University students drop
Data3 <- drop_students(Data3)
table(Data3$University.Student)    # 0
nrow(Data3)                        # 266

fits_B_students <- run_three_lm(Data3)
fit3r1 <- fits_B_students[[1]]; pro5r1 <- fits_B_students[[2]]; pro6r1 <- fits_B_students[[3]]

show_means(Data3)                  # 5.80, 0.71, 0.65

# Resisters (movement dummies = 1) drop
Data4 <- m.data5
nrow(Data4)                        # 318

Data4 <- drop_CA_any1(Data4)
nrow(Data4)                        # 263   

fits_B_resisters <- run_three_lm(Data4)
fit4r1 <- fits_B_resisters[[1]]; pro7r1 <- fits_B_resisters[[2]]; pro8r1 <- fits_B_resisters[[3]]

show_means(Data4)                  # 5.74, 0.71, 0.65

stargazer::stargazer(fit3r1, pro5r1, pro6r1) ## Table 6 of the Supplementary Information
stargazer::stargazer(fit4r1, pro7r1, pro8r1) ## Table 7 of the Supplementary Information


#######################################################################################################################
###### Robustness Checks (2)
#######################################################################################################################
fit.robust1<- lm(Political.Participation.Index ~ Collective.Action+Collective.Action.Violence+Age+Sex+Marriage+
                     Employment+Economic.Status+Location.2+Education+University.Student+
                     +CA1+CA6+CA7+CA8+CA9, data = m.data1, weights = weights); summary(fit.robust1)

fit.robust2 <- lm(Contribution_A~ Collective.Action+Collective.Action.Violence+Age+Sex+Marriage+
                     Employment+Economic.Status+Location.2+Education+University.Student+CA1+CA6+CA7+CA8+CA9, 
                   data = m.data1, weights = weights); summary(fit.robust2)

fit.robust3 <- lm(Contribution_B~ Collective.Action+Collective.Action.Violence+Age+Sex+Marriage+
                     Employment+Economic.Status+Location.2+Education+University.Student+CA1+CA6+CA7+CA8+CA9,  
                   data = m.data1, weights = weights); summary(fit.robust3)

## Table 8 of the Supplementary Information
stargazer(fit.robust1, fit.robust2, fit.robust3) 

#######################################################################################################################
###### Robustness Checks (3) Sensitivity Analyses 
#######################################################################################################################
library(sensemakr)

Sensitivity1 <- sensemakr(model = fit6, 
                          treatment = "Collective.ActionResisters",
                          list(confounder=c("CA6", "CA1")),
                          kd = 1:10,
                          ky = 1:10, 
                          q = 1,
                          alpha = 0.05, 
                          reduce = TRUE)

#### Figure 5 in Supplementary Information 
plot(Sensitivity1)
plot(Sensitivity1, sensitivity.of = "t-value")

Sensitivity2 <- sensemakr(model = fit9, 
                          treatment = "Collective.ActionInstitutional",
                          benchmark_covariates = list(confounder=c("CA6", "CA1")),
                          kd = 1:10,
                          ky = 1:10, 
                          q = 1,
                          alpha = 0.05, 
                          reduce = TRUE)

#### Figure 6 in Supplementary Information 
plot(Sensitivity2)
plot(Sensitivity2, sensitivity.of = "t-value")


#######################################################################################################################
###### Robustness Checks (4) Disaggregating Political Participation Index 
#######################################################################################################################
fit.index1  <- lm(Low.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
                  Employment + Economic.Status + Location.2 + Education + University.Student +
                  CA1 + CA6 + CA7 + CA8 + CA9,
                  data = m.data1, weights = weights)

fit.index2  <- lm(High.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
                  Employment + Economic.Status + Location.2 + Education + University.Student +
                  CA1 + CA6 + CA7 + CA8 + CA9,
                  data = m.data1, weights = weights)

fit.index3  <- lm(Low.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
                  Employment + Economic.Status + Location.2 + Education + University.Student +
                  CA1 + CA6 + CA7 + CA8 + CA9,
                  data = m.data2, weights = weights)

fit.index4  <- lm(High.Participation.Index ~ Collective.Action + Age + Sex + Marriage +
                  Employment + Economic.Status + Location.2 + Education + University.Student +
                  CA1 + CA6 + CA7 + CA8 + CA9,
                  data = m.data2, weights = weights)

## Table 9 in the Supplementary Information 
stargazer(fit.index1, fit.index2, fit.index3, fit.index4)
